Divorce Prediction

lien de la base : https://www.kaggle.com/datasets/vanpatangan/divorce-prediction

1 - Introduction

Le mariage est souvent perçu comme l’union d’une union durable, symbolisant l’engagement et la stabilité dans la vie d’un couple. Pourtant, dans de nombreux contextes, les mariages connaissent des trajectoires variées : certains dure toute une vie, d’autre se terminent plus rapidement par un divorce. Ce phénomène est particulièrement intéressant à observer lorsque celui-ci repose sur un mariage arrangé, qui repose sur des dynamiques sociales et familiales différentes de celles d’un mariage romantique. Ces unions peuvent parfois révéler des différences profondes entre les partenaires ou faire émerger des schémas de relations complexes, voire toxiques.

Dans le cadre de cette étude, nous utilisons une base de données synthétique portant sur des mariages arrangés afin d’examiner la durée de ces unions et les facteurs susceptibles d’influencer leur stabilités. L’Analyse de Survie constitue ici un outil pertinent pour modéliser le temps écoulé entre le mariage et le divorce afin de mieux comprendre la distribution temporelle des ruptures.

Cette analyse est surtout pertinente d’un point de vue sociologique, permet de mieux comprendre les dynamiques relationnelles qui conduisent à la stabilité ou à la rupture d’un mariage. Étudier la durée d’un mariage et les facteurs associés au divorce éclaire notamment des notions essentielles comme la confiance, la communication, la gestion des conflits ou l’évolution des attentes au sein du couple. Comprendre ces mécanismes aide à mieux appréhender la manière dont les individus construisent ou parfois perdent un lien conjugal durable.

Elle présente également un intérêt social plus large : identifier les facteurs de fragilité permet de sensibiliser les couples, d’améliorer les dispositifs d’accompagnement et de renforcer la prévention. L’objectif n’est pas seulement d’anticiper une rupture, mais aussi de favoriser un environnement relationnel plus sain, où les partenaires disposent des ressources pour maintenir un mariage fondé sur la confiance, la solidarité et le respect mutuel.

1.1 Contexte et justification

La stabilité conjugale constitue un enjeu important sur les plans socia, démographique et psychologique. La durée d’un mariage influence notamment le bien-être des individus, le développement des enfants, mais aussi la structuration des familles et la cohésion sociale. À l’inverse, le divorce ou la séparation engendre des conséquences multiples : coûts émotionnels, réorganisation familiale, contraintes économiques ou fragilité psychologique.

Dans le cas des mariages arrangés, ces enjeux sont accentués par des dynamiques culturelles particulières, notamment le rôle de l’entourage, l’absence de choix conjugal initial ou la pression sociale. Étudier la durée de ces unions permet donc de mieux comprendre les mécanismes spécifiques qui favorisent la stabilité ou, au contraire, précipitent la rupture.

1.2 Problématique

Quels facteurs influencent la durée d’un mariage arrangé et la probabilité de divorce ou de séparation au fil du temps ? Comment des caractéristiques individuelles, familiales ou relationnelles peuvent-elles modifier le risque de rupture ?

Quels facteurs influencent la durée de mariage ?

2 Description et préparation des données

2.1 - Presentation des variables

Notre base de données comporte 5000 observations pour 22 variables. Sur les 22 variables, nous retrouvons près de 10 variables quantitatives pour 12 qualitatives. De plus, notre base de données ne comporte aucune valeurs manquantes, ce qui réduit la complexité des prétraitements des données et permet de déterminer directement l’analyse exploratoire. Le tableau ci-dessous synthétise la présentation ainsi que les types et sous-type de variables.

📊 Tableau des variables avec code couleur

Nom_de_la_variable	Description	Type	Sous_type
age_at_marriage	Âge au mariage	Quantitative	Discrète
marriage_duration_years	Durée du mariage	Quantitative	Discrète
divorced	Divorce (oui/non)	Qualitative	Binaire
num_children	Nombre d’enfants	Quantitative	Discrète
education_level	Niveau d’éducation	Qualitative	Ordinale
employment_status	Statut professionnel	Qualitative	Nominale
combined_income	Revenu combiné	Quantitative	Continue
religious_compatibility	Compatibilité religieuse	Qualitative	Nominale
cultural_background_match	Correspondance culturelle	Qualitative	Binaire
communication_score	Score de communication	Quantitative	Continue
conflict_frequency	Fréquence des conflits	Quantitative	Discrète
conflict_resolution_style	Style de résolution de conflit	Qualitative	Nominale
mental_health_issues	Problèmes de santé mentale	Qualitative	Binaire
financial_stress_level	Niveau de stress financier	Quantitative	Continue
infidelity_occurred	Infidélité survenue	Qualitative	Binaire
counseling_attended	A suivi un counseling	Qualitative	Binaire
social_support	Soutien social	Quantitative	Continue
shared_hobbies_count	Nombre de hobbies partagés	Quantitative	Discrète
marriage_type	Type de mariage	Qualitative	Nominale
pre_marital_cohabitation	Cohabitation avant mariage	Qualitative	Binaire
domestic_violence_history	Historique de violence domestique	Qualitative	Binaire
trust_score	Score de confiance	Quantitative	Continue

2.2 - Boxplot des données Quantitatives & Histogrammes

vars_quanti <- Nom_de_la_variable[Type == "Quantitative"]
n <- length(vars_quanti)

# Palette
cols <- brewer.pal(max(n, 3), "Set2")
if (length(cols) < n) cols <- colorRampPalette(cols)(n)

fig <- plot_ly()

for (i in seq_along(vars_quanti)) {

  x <- data.brute[[vars_quanti[i]]]

  # Histogramme
  fig <- fig %>% add_histogram(
    x = x,
    visible = (i == 1),
    xaxis = "x1",
    yaxis = "y1",
    marker = list(color = cols[i]),
    opacity = 0.7,
    showlegend = FALSE
  )

  # Boxplot
  fig <- fig %>% add_boxplot(
    y = x,
    visible = (i == 1),
    xaxis = "x2",
    yaxis = "y2",
    fillcolor = cols[i],
    line = list(color = cols[i]),
    marker = list(color = cols[i]),
    opacity = 0.7,
    boxpoints = "outliers",
    showlegend = FALSE,
    name = ""
  )

}

# Configuration du bouton
buttons <- lapply(seq_len(n), function(i) {

  vis <- rep(FALSE, 2 * n)
  vis[(2*i - 1):(2*i)] <- TRUE

  list(
    method = "update",
    args = list(
      list(visible = vis)
      ),
    label = vars_quanti[i]
  )
})

fig <- fig %>% layout(
  updatemenus = list(
    list(
      type = "dropdown",
      active = 0,
      buttons = buttons,
      x = 0.5,
      y = 1.15,
      xanchor = "center"
    )
  ),
  grid = list(rows = 1, columns = 2, pattern = "independent"),
  xaxis = list(title = "Valeurs"),
  yaxis = list(title = "Effectifs"),
  xaxis2 = list(title = ""),
  yaxis2 = list(title = "Valeurs"),
  width = 700,
  height = 400,
  showlegend = FALSE
)

fig

L’analyse de la variable marriage_duration_years montre une distribution décroissante, avec la majorité des mariages ayant une durée relativement courte. Les effectifs diminuent progressivement lorsque la durée augmente. La durée minimale observée est de 1 ans, la maximale de 40 ans, et la médiane est de 6 ans. On remarque également quelques valeurs extrêmes entre 30 et 40 ans, qui sont isolées par rapport à la majorité des observations. Ces outliers peuvent refléter des cas particuliers de mariages très longs.

2.3 - Histogrammes des données Qualitatives

# Sélection des variables qualitatives
vars_quali <- Nom_de_la_variable[Type == "Qualitative"]
n <- length(vars_quali)

fig <- plot_ly()

for (i in seq_along(vars_quali)) {

  var <- vars_quali[i]
  counts <- table(data.brute[[var]])

  df <- data.frame(
    modalite = names(counts),
    n = as.numeric(counts)
  )

  df$pct <- round(100 * df$n / sum(df$n), 1)

  # Palette pour les modalités
  cols <- brewer.pal(max(3, nrow(df)), "Set2")
  if (length(cols) < nrow(df)) {
    cols <- colorRampPalette(cols)(nrow(df))
  }

  # Construction des PieCharts
  fig <- fig %>% add_pie(
    data = df,
    labels = ~modalite,
    values = ~n,
    visible = (i == 1),
    textinfo = "label+percent",
    marker = list(colors = cols),
    showlegend = FALSE
  )
}

# Configurartion du bouton
buttons <- lapply(seq_len(n), function(i) {

  vis <- rep(FALSE, n)
  vis[i] <- TRUE

  list(
    method = "update",
    args = list(
      list(visible = vis)
    ),
    label = vars_quali[i]
  )
})

## ---- Layout final ----
fig <- fig %>% layout(
  updatemenus = list(
    list(
      type = "dropdown",
      active = 0,
      buttons = buttons,
      x = 0.5,
      y = 1.15,
      xanchor = "center"
    )
  ),
  width = 700,
  height = 500,
  showlegend = FALSE
)

fig

Dans notre dataset, on peut observer un poucentage assez grand de dirvorce en proportion à celui des individus rester mariés, on a 1991 individus sur 5000 qui sont divricés dans notre dataset. Cela représente une proportion de 39,8 % d’individus divorcés. On observe par ailleurs une variable intéressante concernant le type de mariage. Dans notre échantillon, 70,3 % des individus sont mariés par amour, tandis que 24,6 % sont engagés dans un mariage arrangé. Les 5,1 % restants correspondent à d’autres types de mariage, tels que les mariages de convenance, les unions mixtes ou d’autres formes moins fréquentes.

3 - Analyse de survie

Notre base de données comporte une variable temporelle de durée de survie caractérisé par :

• marriage_duration_years : Mesure la Durée du mariage de l’individu.

De plus, nous introduisons une variable \(a\) correspondant à la borne inférieure de la variable de survie. Ici, pour marriage_duration_years, on a \(a = 1\). Cette formalisation permet d’unifier la notation et de clarifier les domaines de définition dans les développements théoriques ultérieurs.

On pose \(X\) la variable aléatoire de survenue de l’évènement d’intérêt, donc le divorce. On note donc les différentes fonctions de survie et leurs interprétations par le tableau suivant :

📊 Interprétation et définitions des fonctions de survie

Fonction	Définition	Durée_du_mariage
\(S(t)\)	\(S(t) = \mathbb{P}(X \gt t) = e^{-H(t)} = e^{-\int_a^t h(u)\,du}\)	Probabilité que le mariage dure ≥ t
\(H(t)\)	\(H(t) = \int_a^t h(u)\,du = -\ln S(t)\)	Risque cumulé de divorce jusqu’à t
\(h(t)\)	\(h(t) = -\dfrac{S'(t)}{S(t)}\)	Risque instantané de divorce à t

Nos données comportent une censure : certains individus n’ont pas encore connu l’événement d’intérêt, c’est à dire qu’ils sont toujours encore mariés. Cette information est déjà inscrite dans la base de données via la variable divorced, qui indique si l’individu est divorcé ou non que l’on note :

\[ \delta_i = \begin{cases} 1 & \text{si l'événement divorce est observé pour } i \\ 0 & \text{si l'observation n'est pas divorcé} \end{cases} \]

3.1 Censure

Soit \(X_i\) le temps de survie réel de l’individu \(i\) (durée jusqu’à l’événement d’intérêt, ici le divorce), et \(C_i\) la variable aléatoire du temps de censure, représentant le moment auquel l’individu quitte l’étude ou n’a pas encore subi l’événement.

La durée réellement observée pour chaque individu dépend du type de censure :

• Censure à droite : \[ T_i = \min(X_i, C_i) \]
• Censure à gauche : \[ T_i = \max(X_i, C_i) \]

3.1.1 Censure à droite

La censure à droite se produit lorsqu’un individu n’a pas encore subi l’événement d’intérêt (ici le divorce) au moment de sa dernière observation (\(X_i > C_i\)).
Les principaux types de censure à droite sont :

• Type I (fixe) : tous les individus sont censurés au même moment prédéterminé.
  
• Type II : l’étude s’arrête dès qu’un certain nombre d’événements est observé.
  
• Type III (aléatoire) : le moment de censure varie d’un individu à l’autre, par exemple à cause de la fin de suivi variable, de pertes de vue ou d’arrêts de participation. Ce type est le plus courant dans les études observationnelles.

3.1.2 Censure à gauche

La censure à gauche se produit lorsque l’événement a eu lieu avant le début de l’observation, et on ne connaît que la borne supérieure du temps de survie (\(X_i < C_i\)).
Elle est beaucoup plus rare dans les études humaines et moins souvent traitée dans la littérature.

3.1.3 Censure par intervalle

Une censure par intervalle survient lorsqu’on sait seulement que l’événement s’est produit entre deux dates d’observation. Dans la pratique, elle est souvent convertie en censure à droite pour simplifier l’analyse.

---

Dans notre base de données, certains mariages n’ont pas abouti à un divorce au moment de la fin de l’étude, et le temps de suivi varie selon les individus.
On en déduit que les données présentent une censure à droite de type III (aléatoire).
On suppose que cette censure est non informative, c’est-à-dire indépendante de la probabilité de divorce, conformément aux hypothèses classiques des modèles de survie.

Dans ce contexte, la durée réellement observée pour chaque mariage est donnée par :

\[ T = \min(X, C) \]

3.2 - Estimateur de la Fonction de survie S(t)

Estimateur empirique de la fonction de survie :

\[ \hat{S}(t) = \frac{1}{n} \sum_{i=1}^{n} \boldsymbol{1}_{\{t_i > t\}} \]

• \(n\) = nombre total d’observations
  
• \(t_i\) = temps jusqu’à l’événement pour l’individu \(i\)
  
• \(\boldsymbol{1}_{\{t_i > t\}}\) = indicateur qui vaut 1 si l’individu n’a pas encore eu l’événement à \(t\), 0 sinon

Cet estimateur correspond simplement à la proportion d’individus encore mariés au temps \(t\).
Il suppose qu’il n’y a aucune donnée censurée, c’est-à-dire que tous les individus ont eu l’événement observé.

📊 Estimateur de survie sans censure, variance binomiale et intervalle de confiance plain à 95 %

Méthode	Formule	Description
Estimateur empirique de survie (sans censure)	\(\hat{S}(t)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{\{t_i\gt t\}}\)	Dans le cas sans censure, Kaplan–Meier coïncide avec l’estimateur empirique de la fonction de survie.
Variance (loi binomiale, cas sans censure)	\(\widehat{\text{Var}}[\hat{S}(t)] = \frac{\hat S(t) (1 - \hat S(t))}{n}\)	Variance estimée selon la loi binomiale, adaptée aux données entièrement observées.
Intervalle de confiance plain à 95 %	\(\text{IC}_{95\%}(t) = \hat S(t) \pm 1.96 \sqrt{\widehat{\text{Var}}[\hat S(t)]}\)	Intervalle de confiance classique basé sur la variance binomiale.

---

3.3 - Estimateur avec censure (Kaplan-Meier)

L’estimateur de Kaplan-Meier découle de l’idée suivante : survivre après un temps \(t_n\) revient à être vivant juste avant \(t_n\) et ne pas subir l’événement à ce temps. Formellement, pour \(t_0 < t_1 < \dots < t_{n-1} < t_n\) :

La probabilité de survie jusqu’à \(t_n\) peut s’écrire en utilisant la règle de multiplication des probabilités :

\[ \mathbb{P}(X > t_n) = \mathbb{P}(X > t_1, X > t_2, \dots, X > t_n) \]

On introduit une récurrence : pour tout \(k \ge 1\),

\[ \mathbb{P}(X > t_k \mid X > t_{k-1}, \dots, X > t_1) = \mathbb{P}(X > t_k \mid X > t_{k-1}) \]

où l’égalité découle de l’indépendance conditionnelle induite par l’ordre croissant des temps.

Ainsi, par récurrence sur les indices \(k\) :

\[ \begin{aligned} \mathbb{P}(X > t_1, X > t_2, \dots, X > t_n) &= \mathbb{P}(X > t_1) \cdot \mathbb{P}(X > t_2 \mid X > t_1) \\ &\quad \cdot \mathbb{P}(X > t_3 \mid X > t_1, X > t_2) \cdots \mathbb{P}(X > t_n \mid X > t_1, \dots, X > t_{n-1}) \\ &= \mathbb{P}(X > t_1) \prod_{k=2}^{n} \mathbb{P}(X > t_k \mid X > t_{k-1}) \end{aligned} \]

On considère les temps d’événements distincts \(T_{(1)} < T_{(2)} < \dots < T_{(j)}\) (décès ou divorce observés) rangés par ordre croissant.
On définit \(T_{(0)} = 0\), la borne inférieure du temps (par exemple \(a=1\) pour la durée de mariage).

Ainsi, la probabilité de survie jusqu’au temps \(T_{(j)}\) peut s’écrire comme un produit de probabilités conditionnelles :

\[ \begin{aligned} \mathbb{P}(X > T_{(j)}) &= \prod_{k=1}^{j} \mathbb{P}(X > T_{(k)} \mid X > T_{(k-1)}) \end{aligned} \]

Pour chaque temps d’événement \(T_{(k)}\), on s’intéresse à la probabilité conditionnelle de subir l’événement à ce temps, sachant que l’individu était encore à risque juste avant :

\[ \mathbb{P}(X \le T_{(k)} \mid X > T_{(k-1)}) \]

Cette quantité représente la probabilité qu’un individu qui a « survécu » jusqu’à \(T_{(k-1)}\) subisse l’événement à \(T_{(k)}\).

En pratique, on dispose des données observées :

• \(n_k\) = nombre d’individus encore à risque juste avant \(T_{(k)}\)
  
• \(d_k\) = nombre d’événements observés à \(T_{(k)}\)

On peut alors estimer cette probabilité conditionnelle par :

\[ \hat{\mathbb{P}}(X \le T_{(k)} \mid X > T_{(k-1)}) = \frac{d_k}{n_k} \]

La probabilité de survivre au temps \(T_{(k)}\) est le complémentaire :

\[ \hat{q}_k = \hat{\mathbb{P}}(X \ge T_{(k)} \mid X > T_{(k-1)}) = 1 - \hat{\mathbb{P}}(X \le T_{(k)} \mid X > T_{(k-1)}) = 1 - \frac{d_k}{n_k} \]

Enfin, en remplaçant les probabilités conditionnelles dans le produit de survie, on obtient l’estimateur de Kaplan-Meier (ou produit-limite) :

\[ \hat{S}(t) = \prod_{T_{(k)} \le t} \hat{q}_k = \prod_{T_{(k)} \le t} \left( 1 - \frac{d_k}{n_k} \right) \]

Ainsi, l’estimateur de Kaplan-Meier corrige naturellement le biais dû à la censure et fournit une estimation non paramétrique de la fonction de survie.

📊 Estimateur de Kaplan-Meier, variance de Greenwood et intervalle de confiance log à 95 %

Méthode	Formule	Description
Kaplan-Meier	\(\hat{S}(t) = \prod_{T_{(k)} \le t} \left( 1 - \dfrac{d_k}{n_k} \right)\)	Estimateur non paramétrique de la fonction de survie basé sur les événements observés et le nombre d’individus à risque.
Variance de Greenwood	\(\widehat{\text{Var}}\left[\hat{S}(t)\right] = \hat{S}(t)^2 \sum_{T_{(k)} \le t} \dfrac{d_k}{n_k (n_k - d_k)}\)	Variance estimée de Kaplan-Meier selon la formule de Greenwood.
Intervalle de confiance log à 95 %	\(\text{IC}_{95\%}(t) = \hat S(t) \pm 1.96 \sqrt{\widehat{\text{Var}}[\hat S(t)]}\)	Intervalle de confiance construit via une transformation logarithmique de S(t), qui est la méthode ‘plain’ de survfit().

3.3.1 - Graphique Kaplan-Meier

# Données
time <- data.brute$marriage_duration_years
event <- data.brute$divorced

# Survie sans censure avec IC
surv_no_censor <- survfit(Surv(time, rep(1, length(time))) ~ 1, conf.type = "plain")
surv_no_censor$lower[is.na(surv_no_censor$lower)] <- 0
surv_no_censor$upper[is.na(surv_no_censor$upper)] <- 0

n_no <- length(surv_no_censor$time)

# Survie avec censure avec IC
surv_censor <- survfit(Surv(time, event) ~ 1, conf.type = "plain")

# Plotly avec IC et legendgroup
fig <- plot_ly() %>%
  
  # Sans censure
  add_ribbons(x = surv_no_censor$time,
              ymin = surv_no_censor$lower,
              ymax = surv_no_censor$upper,
              fillcolor = "#1f78b433",
              line = list(color = "transparent"),
              legendgroup = "nocens",
              showlegend = FALSE) %>%
  add_lines(x = surv_no_censor$time,
            y = surv_no_censor$surv,
            name = "Sans censure",
            line = list(color = "#1f78b4", width = 2),
            legendgroup = "nocens",
            showlegend = TRUE) %>%
  
  # Avec censure
  add_ribbons(x = surv_censor$time,
              ymin = surv_censor$lower,
              ymax = surv_censor$upper,
              fillcolor = "#e31a1c33",
              line = list(color = "transparent"),
              legendgroup = "cens",
              showlegend = FALSE) %>%
  add_lines(x = surv_censor$time,
            y = surv_censor$surv,
            name = "Avec censure",
            line = list(color = "#e31a1c", width = 2),
            legendgroup = "cens",
            showlegend = TRUE) %>%
  
  # Layout
  layout(title = "Comparaison des fonctions de survie : sans et avec censure",
         xaxis = list(title = "Durée du mariage (années)"),
         yaxis = list(title = "Probabilité de rester marié"),
         legend = list(x = 0.80, y = 0.95))

div(style="text-align:center;", fig)

Pour la courbe sans censure, la probabilité de rester marié diminue progressivement avec l’augmentation de la durée du mariage. Par exemple, après 1 an de mariage, les couples ont environ 82 % de rester mariés, mais cette proportion tombe à 32 % après 10 ans et approche de zéro après 36 à 40 ans. Les intervalles de confiance sont étroits au début car beaucoup de couples sont encore à risque. Comme les individus censurés ne sont pas pris en compte, cette estimation sous‑estime la survie réelle.

Pour la courbe avec censure, la probabilité de rester marié est plus élevée car les couples censurés (par exemple ceux pour lesquels on ne connaît pas la fin du mariage) ne sont pas considérés comme ayant divorcé. Après 1 an de mariage, les couples ont une chance de 93 % de rester mariés. Cette proportion descend à 64 % après 10 ans et à environ 10 % après 40 ans. Les intervalles de confiance s’élargissent légèrement avec la durée du mariage car moins de couples restent à risque. Cette estimation correspond à l’estimateur de Kaplan‑Meier et reflète mieux la survie réelle des mariages dans la population étudiée.

En résumé, sans censure, la courbe montre la survie brute et sous-estime la durée réelle des mariages, tandis que la courbe avec censure ajuste pour les mariages dont la fin n’a pas été observée, donnant une estimation plus fiable de la probabilité de rester marié dans le temps.

3.4 Estimateur de survie avec covariables

Dans cette section, nous évaluons la fonction de survie en fonction de différentes variables explicatives afin de déterminer si certains groupes ont une influence sur la courbe de survie de Kaplan–Meier.

Nous réaliserons des tests d’hypothèses pour vérifier si les courbes de survie diffèrent significativement entre les groupes.

En particulier, nous utiliserons des tests tels que le Log-Rank ou le test de Petro & Prentice, selon le type de covariable étudiée.

3.4.1 Cas général de comparaison de plusieurs groupes

On considère \(k\) groupes de survie :

\[ S_1(t),\dots,S_k(t) \]

Hypothèses globales :

\[ \begin{cases} H_0 : S_1(t)=\dots=S_k(t), & \forall t \\ H_1 : \exists r,s,t \text{ tels que } S_r(t)\neq S_s(t) \end{cases} \]

Soient les temps distincts de décès :

\[ T_1 < \dots < T_N \]

Pour chaque temps \(T_i\) et chaque groupe \(g=1,\dots,k\) :

• \(d_{gi}\) = nombre de décès observés dans le groupe \(g\) au temps \(T_i\)
• \(n_{gi}\) = nombre d’individus à risque dans le groupe \(g\) juste avant \(T_i\)

Sommes sur les groupes :

\[ d_i = \sum_{g=1}^k d_{gi}, \quad n_i = \sum_{g=1}^k n_{gi} \]

Variables aléatoires associées :

\[ D_{gi} \text{ dont la valeur observée est } d_{gi} \]

On empile les \(k\) nombres de décès observés au temps \(i\) :

\[ \mathbf{D}_i = \begin{pmatrix} D_{1i}\\ \vdots\\ D_{ki} \end{pmatrix} \in \mathbb{R}^k \]

Sous \(H_0\), l’espérance conditionnelle :

\[ \mathbb{E}(\mathbf{D}_i) = \frac{d_i}{n_i} \begin{pmatrix} n_{1i} \\ \vdots \\ n_{ki} \end{pmatrix}, \quad \mathbb{E}(D_{gi}) = \frac{n_{gi} d_i}{n_i} \]

\[ \mathbf{V}_i = \mathbb{V}(\mathbf{D}_i) = \frac{n_i - d_i}{n_i - 1} \cdot \frac{d_i}{n_i^2} \Big( \begin{pmatrix} n_{1i} & n_{2i} & \dots & n_{ki} \end{pmatrix} I_k - \frac{1}{n_i} \begin{pmatrix} n_{1i} & n_{2i} & \dots & n_{ki} \end{pmatrix}^\top \begin{pmatrix} n_{1i} & n_{2i} & \dots & n_{ki} \end{pmatrix} \Big) \]

Vecteur score log-rank généralisé :

\[ \mathbf{U} = \sum_{i=1}^N w_i (\mathbf{D}_i - \mathbb{E}(\mathbf{D}_i)) \in \mathbb{R}^k \]

Matrice de variance :

\[ \mathbf{V} = \sum_{i=1}^N w_i^2 \mathbf{V}_i \in \mathbb{R}^{k\times k} \]

Statistique de test :

\[ \chi^2 = \mathbf{U}^\top \mathbf{V}^{-1} \mathbf{U} \sim \chi^2_{k-1} \]

Ainsi selon le test, on a :

• Log-Rank lorsque pour tout i \(\in\) [1, N], \(w_i = 1\)
  
• Petro & Prentice lorsque \(w_i = \prod_{k=1}^{i} \frac{n_k}{n_k + d_k}\)

3.4.2 Exemple selon deux groupes (k = 2)

# Données
time1 <- data.brute$marriage_duration_years[data.brute$mental_health_issues==0]
event1 <- data.brute$divorced[data.brute$mental_health_issues==0]

time2 <- data.brute$marriage_duration_years[data.brute$mental_health_issues==1]
event2 <- data.brute$divorced[data.brute$mental_health_issues==1]


# Survie avec censure avec IC
surv_censor1 <- survfit(Surv(time1, event1) ~ 1, conf.type = "plain")
surv_censor2 <- survfit(Surv(time2, event2) ~ 1, conf.type = "plain")

# Plotly avec IC et legendgroup
fig <- plot_ly() %>%

  # sans infidéité
  add_ribbons(x = surv_censor1$time,
              ymin = surv_censor1$lower,
              ymax = surv_censor1$upper,
              fillcolor = "#e31a1c33",
              line = list(color = "transparent"),
              legendgroup = "noinfi",
              showlegend = FALSE) %>%
  add_lines(x = surv_censor1$time,
            y = surv_censor1$surv,
            name = "Sans problème de santé mentale",
            line = list(color = "#e31a1c", width = 2),
            legendgroup = "noinfi",
            showlegend = TRUE) %>%
  
  # Avec infidélité
  add_ribbons(x = surv_censor2$time,
              ymin = surv_censor2$lower,
              ymax = surv_censor2$upper,
              fillcolor = "#1f78b433",
              line = list(color = "transparent"),
              legendgroup = "infi",
              showlegend = FALSE) %>%
  add_lines(x = surv_censor2$time,
            y = surv_censor2$surv,
            name = "Avec problème de santé mentale",
            line = list(color = "#1f78b4", width = 2),
            legendgroup = "infi",
            showlegend = TRUE) %>%
  
  # --- Layout ---
  layout(title = "Fonctions de survie de Kaplan Meier sans et avec problème de santé mentale",
         xaxis = list(title = "Durée du mariage (années)"),
         yaxis = list(title = "Probabilité de rester marié"),
         legend = list(x = 0.80, y = 0.95))

div(style="text-align:center;", fig)

Dans notre cas de figure, au dessus, on pose deux groupes :

• Groupe 1 : sans problème de santé mentale
  
• Groupe 2 : avec problème de santé mentale

On évoque donc les hypotèses suivantes :

\[ \begin{cases} H_0 : S_1(t) = S_2(t), & \forall t \\ H_1 : S_1(t) \neq S_2(t), & \exists t \end{cases} \]

Soit i \(\in\) [1, 40], on a :

Vecteur de divorce :

\[ \mathbf{D}_i = \begin{pmatrix} D_{1i} \\ D_{2i} \end{pmatrix}, \quad \mathbb{E}(\mathbf{D}_i) = \frac{d_i}{n_i} \begin{pmatrix} n_{1i} \\ n_{2i} \end{pmatrix} \]

Variance du premier composant (groupe 1) :

\[ \mathbb{V}(D_{1i}) = \frac{(n_i - d_i)}{n_i - 1} \frac{d_i n_{1i} n_{2i}}{n_i^2} \]

Vecteur score réduit à une dimension :

\[ U = \sum_{i=1}^{40} w_i (D_{1i} - E(D_{1i})) \]

Variance :

\[ \text{Var}(U) = \sum_{i=1}^{40} w_i^2 \mathbb{V}(D_{1i}) \]

Statistique de test log-rank :

\[ \chi_0^2 = \frac{U^2}{\text{Var}(U)} \sim \chi_1^2 \]

On effectue le test selon le test de log-rank :

# vecteurs concaténés
time  <- c(time1, time2)
event <- c(event1, event2)

# facteur groupe
group <- factor(c(
  rep("no_mental_health_issue", length(time1)),
  rep("mental_health_issue", length(time2))
))

# test du log-rank (khi-deux)
res <- survdiff(Surv(time, event) ~ group, rho = 0)
res

## Call:
## survdiff(formula = Surv(time, event) ~ group, rho = 0)
## 
##                                 N Observed Expected (O-E)^2/E (O-E)^2/V
## group=mental_health_issue    1019      443      388      7.70      10.1
## group=no_mental_health_issue 3981     1548     1603      1.86      10.1
## 
##  Chisq= 10.1  on 1 degrees of freedom, p= 0.002

On effectue le test selon le test de Peto & Prentice :

# vecteurs concaténés
time  <- c(time1, time2)
event <- c(event1, event2)

# facteur groupe
group <- factor(c(
  rep("no_mental_health_issue", length(time1)),
  rep("mental_health_issue", length(time2))
))

# test du log-rank (khi-deux)
res <- survdiff(Surv(time, event) ~ group, rho = 1)
res

## Call:
## survdiff(formula = Surv(time, event) ~ group, rho = 1)
## 
##                                 N Observed Expected (O-E)^2/E (O-E)^2/V
## group=mental_health_issue    1019      328      294      4.08      6.62
## group=no_mental_health_issue 3981     1143     1178      1.02      6.62
## 
##  Chisq= 6.6  on 1 degrees of freedom, p= 0.01

Le test du Log-Rank (rho = 0) donne une statistique de Chi-2 égale à 10,1 avec 1 degré de liberté et une p-value de 0,002. Comme la p-value est inférieure au seuil de 0,05, l’hypothèse nulle \(H_0\) d’égalité des fonctions de survie entre les deux groupes peut être rejetée. On en conclut qu’il existe une différence significative entre la survie des individus avec et sans problème de santé mentale.

Le test de Petro & Prentice (rho = 1), qui pondère davantage les événements précoces, fournit une statistique de Chi-2 de 6,62 avec 1 degré de liberté et une p-value de 0,01. Cette p-value étant également inférieure à 0,05, le résultat confirme que la différence observée entre les groupes est significative, même en mettant un poids plus important sur les événements survenus tôt.

3.4.3 Exemple selon plusieurs groupes (k = 4 > 2)

# Données
time1 <- data.brute$marriage_duration_years[data.brute$education_level=="Bachelor"]
event1 <- data.brute$divorced[data.brute$education_level=="Bachelor"]

time2 <- data.brute$marriage_duration_years[data.brute$education_level=="High School"]
event2 <- data.brute$divorced[data.brute$education_level=="High School"]

time3 <- data.brute$marriage_duration_years[data.brute$education_level=="Master"]
event3 <- data.brute$divorced[data.brute$education_level=="Master"]

time4 <- data.brute$marriage_duration_years[data.brute$education_level=="PhD"]
event4 <- data.brute$divorced[data.brute$education_level=="PhD"]

time5 <- data.brute$marriage_duration_years[data.brute$education_level=="No Formal Education"]
event5 <- data.brute$divorced[data.brute$education_level=="No Formal Education"]


# Survie avec censure avec IC
surv_censor1 <- survfit(Surv(time1, event1) ~ 1, conf.type = "plain")
surv_censor2 <- survfit(Surv(time2, event2) ~ 1, conf.type = "plain")
surv_censor3 <- survfit(Surv(time3, event3) ~ 1, conf.type = "plain")
surv_censor4 <- survfit(Surv(time4, event4) ~ 1, conf.type = "plain")
surv_censor5 <- survfit(Surv(time5, event5) ~ 1, conf.type = "plain")

cols <- c(
  "Bachelor" = "#1f77b4",
  "High School" = "#e31a1c",
  "Master" = "#33a02c",
  "PhD" = "#ff7f00",
  "No Formal Education" = "#6a3d9a"
)


# Plotly avec IC et legendgroup
fig <- plot_ly() %>%
  
  # Bachelor
  add_ribbons(
    x = surv_censor1$time,
    ymin = surv_censor1$lower,
    ymax = surv_censor1$upper,
    fillcolor = paste0(cols["Bachelor"], "33"),
    line = list(color = "transparent"),
    legendgroup = "Bachelor",
    showlegend = FALSE
  ) %>%
  add_lines(
    x = surv_censor1$time,
    y = surv_censor1$surv,
    name = "Bachelor",
    line = list(color = cols["Bachelor"], width = 2),
    legendgroup = "Bachelor",
    showlegend = TRUE
  ) %>%
  
  # High School
  add_ribbons(
    x = surv_censor2$time,
    ymin = surv_censor2$lower,
    ymax = surv_censor2$upper,
    fillcolor = paste0(cols["High School"], "33"),
    line = list(color = "transparent"),
    legendgroup = "High School",
    showlegend = FALSE
  ) %>%
  add_lines(
    x = surv_censor2$time,
    y = surv_censor2$surv,
    name = "High School",
    line = list(color = cols["High School"], width = 2),
    legendgroup = "High School",
    showlegend = TRUE
  ) %>%
  
  # Master
  add_ribbons(
    x = surv_censor3$time,
    ymin = surv_censor3$lower,
    ymax = surv_censor3$upper,
    fillcolor = paste0(cols["Master"], "33"),
    line = list(color = "transparent"),
    legendgroup = "Master",
    showlegend = FALSE
  ) %>%
  add_lines(
    x = surv_censor3$time,
    y = surv_censor3$surv,
    name = "Master",
    line = list(color = cols["Master"], width = 2),
    legendgroup = "Master",
    showlegend = TRUE
  ) %>%
  
  # PhD
  add_ribbons(
    x = surv_censor4$time,
    ymin = surv_censor4$lower,
    ymax = surv_censor4$upper,
    fillcolor = paste0(cols["PhD"], "33"),
    line = list(color = "transparent"),
    legendgroup = "PhD",
    showlegend = FALSE
  ) %>%
  add_lines(
    x = surv_censor4$time,
    y = surv_censor4$surv,
    name = "PhD",
    line = list(color = cols["PhD"], width = 2),
    legendgroup = "PhD",
    showlegend = TRUE
  ) %>%
  
  # No formal education
  add_ribbons(
    x = surv_censor5$time,
    ymin = surv_censor5$lower,
    ymax = surv_censor5$upper,
    fillcolor = paste0(cols["No Formal Education"], "33"),
    line = list(color = "transparent"),
    legendgroup = "No Formal Education",
    showlegend = FALSE
  ) %>%
  add_lines(
    x = surv_censor5$time,
    y = surv_censor5$surv,
    name = "No formal education",
    line = list(color = cols["No Formal Education"], width = 2),
    legendgroup = "No Formal Education",
    showlegend = TRUE
  ) %>%
  
  layout(
    title = "Fonctions de survie par niveau d'éducation",
    xaxis = list(title = "Durée du mariage (années)"),
    yaxis = list(title = "Probabilité de rester marié"),
    legend = list(x = 0.80, y = 0.95)
  )

div(style="text-align:center;", fig)

Dans notre cas de figure, au dessus, on pose deux groupes :

• Groupe 1 : Bachelor
• Groupe 2 : High School
• Groupe 3 : Master
• Groupe 4 : PhD

On évoque donc les hypotèses suivantes :

\[ \begin{cases} H_0 : S_1(t) = S_2(t) = S_3(t) = S_4(t) ,\ \forall t \\ H_1 : \exists t,\ \exists (i,j)\in\{1,\dots,4\}^2,\ i \neq j,\ \text{tel que } S_i(t)\neq S_j(t) \end{cases} \]

On obtient donc selon le test du log-rank :

# vecteurs concaténés
time  <- c(time1, time2, time3, time4)
event <- c(event1, event2, event3, event4)

# facteur groupe
group <- factor(c(
  rep("Bachelor", length(time1)),
  rep("Hight School", length(time2)),
  rep("Master", length(time3)),
  rep("PhD", length(time4))
))

# test du log-rank (khi-deux)
res <- survdiff(Surv(time, event) ~ group, rho = 0)
res

## Call:
## survdiff(formula = Surv(time, event) ~ group, rho = 0)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## group=Bachelor     2069      802    833.8    1.2139     2.271
## group=Hight School 1513      605    591.5    0.3087     0.471
## group=Master        963      402    396.4    0.0786     0.105
## group=PhD           224      102     89.3    1.8127     2.013
## 
##  Chisq= 3.6  on 3 degrees of freedom, p= 0.3

On obtient donc selon le test de Peto & Pretice :

# vecteurs concaténés
time  <- c(time1, time2, time3, time4)
event <- c(event1, event2, event3, event4)

# facteur groupe
group <- factor(c(
  rep("Bachelor", length(time1)),
  rep("Hight School", length(time2)),
  rep("Master", length(time3)),
  rep("PhD", length(time4))
))

# test du log-rank (khi-deux)
res <- survdiff(Surv(time, event) ~ group, rho = 1)
res

## Call:
## survdiff(formula = Surv(time, event) ~ group, rho = 1)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## group=Bachelor     2069    593.6    616.2   0.82612   1.92530
## group=Hight School 1513    450.0    440.1   0.22186   0.42204
## group=Master        963    288.0    288.7   0.00172   0.00285
## group=PhD           224     78.3     64.9   2.75899   3.80173
## 
##  Chisq= 5  on 3 degrees of freedom, p= 0.2

Le test du Log-Rank (rho = 0) donne une statistique de Chi-2 égale à 3,6 avec 3 degrés de liberté et une p-value de 0,3. Comme la p-value est supérieure au seuil de 0,05, l’hypothèse nulle H0 d’égalité des fonctions de survie entre les niveaux d’éducation ne peut pas être rejetée. On en conclut qu’il n’existe pas de différence significative de survie selon le niveau d’éducation.

Le test de Peto & Prentice (rho = 1), qui pondère davantage les événements précoces, fournit une statistique de Chi-2 de 5 avec 3 degrés de liberté et une p-value de 0,2. Cette p-value étant également supérieure à 0,05, le résultat confirme qu’aucune différence significative de survie n’est observée entre les niveaux d’éducation, même en mettant un poids plus important sur les événements survenus tôt.

3.5 - Estimateur de Nelson-Aalen

L’estimateur de Nelson-Aalen permet d’estimer le risque cumulatif \(h(t)\) dans le cadre de données censurées.

On définit :

• \(H(t) = \mathbb{P}(T > t) = \mathbb{P}(X > t, C > t) = \mathbb{P}(X > t)\mathbb{P}(C > t)= S(t) G(t)\) où \(G\) est la fonction de survie de la censure \(C\)

• \(H_1(t) = \mathbb{P}(T > t, \delta = 1) = \mathbb{P}(X > t, C > X)\)

On peut écrire \(H_1(t)\) en fonction de la densité \(f(u)\) de \(X\) et de \(G(u)\) :

\[ \begin{aligned} H_1(t) &= \mathbb{P}(X > t,\, C > X) \\ &= \mathbb{E}\big[ \mathbf{1}_{\{X > t\}} \cdot \mathbf{1}_{\{C > X\}} \big] \\[6pt] &= \mathbb{E}\Big[ \mathbf{1}_{\{X > t\}} \, \mathbb{E}\big[\mathbf{1}_{\{C > X\}}\mid X\big] \Big] \\[6pt] &= \mathbb{E}\big[ \mathbf{1}_{\{X > t\}} \, \mathbb{P}(C > X \mid X) \big] \\[6pt] &= \mathbb{E}\big[ \mathbf{1}_{\{X > t\}} \, G(X^-) \big] \\[6pt] &= \displaystyle \int_{t}^{\infty} G(u^-) \, f(u)\,du \\[6pt] &= - \displaystyle \int_{t}^{\infty} G(u^-) \, dS(u) \end{aligned} \]

On obtient donc :

\[ dH_1(t) = G(t^{-})dS(t) \]

Et donc par le temps on obtient :

\[ \frac{dH_1(t)}{dt} = \frac{G(t^{-})dS(t)}{dt} \]

ce qui donne mathématiquement :

\[ H_1'(t) = G(t^{-})S'(t) \]

Ainsi on a :

\[ \begin{aligned} \hat{H}_{NA}(t) &= \displaystyle \int_{0}^{t} h(u) \, du \\[2mm] &= \displaystyle \int_{0}^{t} -\frac{S'(u)}{S(u)} \, du \\[2mm] &= \displaystyle \int_{0}^{t} -\frac{\frac{H_1(u)}{G(u^{-})}}{\frac{H(u)}{G(u)}} \, du \\[2mm] &= \displaystyle \int_{0}^{t} -\frac{H_1(u)}{H(u)}\frac{G(u)}{G(u^{-})} \, du \\[2mm] &= \displaystyle \int_{0}^{t} -\frac{H_1(u)}{H(u)} \, du \end{aligned} \]

Un estimateur naturel s’obtient en remplaçant les fonctions \(H\) et \(H_1\) par leurs équivalents empiriques (calculables car les variables \(T\) et \(\delta\) sont observées):

\[ \hat{H}(u) = \frac{1}{n} \sum_{i=1}^{n} \mathbf{1}_{\{T_i > u\}}, \quad \hat{H}_1(u) = \frac{1}{n} \sum_{i=1}^{n} \mathbf{1}_{\{T_i > u, \delta_i = 1\}} \]

L’estimateur de Nelson-Aalen est alors donné par :

\[ \hat{H}_{NA}(t) = \displaystyle \int_{0}^{t} - \frac{\displaystyle \sum_{i=1}^{n} \mathbf{1}_{\{T_i > u, \delta_i = 1\}}}{\displaystyle \sum_{i=1}^{n} \mathbf{1}_{\{T_i > u\}}} \, du \]

Comme \(T\) est à temps discret, l’intégrale devient une somme sur les temps d’événement distincts , et on définit alors pour chaque temps d’événement \(t_i\) :

\[ d_i = \sum_{j=1}^{n} \mathbf{1}_{\{T_j = t_i, \delta_j = 1\}}, \quad n_i = \sum_{j=1}^{n} \mathbf{1}_{\{T_j \ge t_i\}}. \]

Ce qui donne :

\[ \hat{H}_{NA}(t) = \sum_{t_i \le t} \frac{d_i}{n_i}. \]

---

Une autre façon de calculer la fonction de risque cumulée et de passer par l’estimateur de beslow.

Rappel : l’estimateur de Kaplan–Meier de la fonction de survie s’écrit, pour des temps d’événement distincts \(t_1<\dots<t_m\), \[ \hat{S}(t)=\prod_{t_i\le t}\left(1-\frac{d_i}{n_i}\right), \] où \(d_i\) est le nombre d’événements au temps \(t_i\) et \(n_i\) le nombre d’individus à risque juste avant \(t_i\).

En utilisant la relation \[ H(t)=-\log S(t), \] on obtient l’estimateur de Breslow du risque cumulé : \[ \hat{H}_{\text{Breslow}}(t) = -\log\big(\hat{S}(t)\big) = -\log\!\left(\prod_{t_i\le t}\left(1-\frac{d_i}{n_i}\right)\right) = -\sum_{t_i\le t} \log\!\left(1-\frac{d_i}{n_i}\right). \]

Pour des fractions \(d_i/n_i\) petites, on utilise l’approximation \(\log(1-x)\approx -x\) pour \(x\) proche de \(0\). Ainsi \[ \sum_{t_i\le t}\log\!\left(1-\frac{d_i}{n_i}\right) \approx \sum_{t_i\le t}\frac{d_i}{n_i}, \] Ce qui montre que l’estimateur de Breslow est proche (et asymptotiquement équivalent) à l’estimateur de Nelson–Aalen \(\hat{H}_{NA}(t)=\sum_{t_i\le t}\dfrac{d_i}{n_i}\) lorsque les sauts sont petits.

📊 Méthodes d’estimation du risque cumulatif : Nelson-Aalen et Breslow

Méthode	Formule	Variance	Description
Nelson-Aalen	\(\hat{H}_{NA}(t) = \sum_{t_k \le t} \dfrac{d_k}{n_k}\)	\(\text{Var}(\hat{H}_{NA}(t)) = \sum_{t_k \le t} \dfrac{d_k}{n_k^2}\)	Estimateur non paramétrique basé sur les événements observés et le nombre de sujets à risque.
Breslow	\(\hat{H}_{\text{Breslow}}(t) = - \sum_{t_k \le t} \log\left(1 - \dfrac{d_k}{n_k}\right)\)	\(\text{Var}(\hat{H}_{\text{Breslow}}(t)) = \sum_{t_k \le t} \dfrac{d_k}{n_k(n_k - d_k)}\)	Estimateur du risque cumulatif dérivé de \(H(t) = -\log(S(t))\) via l’estimateur de Kaplan-Meier.

3.5.1 Graphique de la fonction du risque cumulatif

# Données 
time  <- data.brute$marriage_duration_years
event <- data.brute$divorced

# Avec censure 
km_cens <- survfit(Surv(time, event) ~ 1)

# Nelson-Aalen
H_NA_cens <- cumsum(km_cens$n.event / km_cens$n.risk)
se_H_NA_cens <- sqrt(cumsum(km_cens$n.event / (km_cens$n.risk^2)))
df_NA_cens <- data.frame(
  time = km_cens$time,
  H = H_NA_cens,
  lower = H_NA_cens - 1.96 * se_H_NA_cens,
  upper = H_NA_cens + 1.96 * se_H_NA_cens
)

# Breslow
km_bres_cens <- survfit(Surv(time, event) ~ 1, type = "kaplan-meier")
df_Bres_cens <- data.frame(
  time = km_bres_cens$time,
  H = -log(km_bres_cens$surv),
  lower = -log(km_bres_cens$upper),
  upper = -log(km_bres_cens$lower)
)

# Sans censure 
km_nocens <- survfit(Surv(time, rep(1,length(time))) ~ 1)

# Nelson-Aalen
H_NA_nocens <- cumsum(km_nocens$n.event / km_nocens$n.risk)
se_H_NA_nocens <- sqrt(cumsum(km_nocens$n.event / (km_nocens$n.risk^2)))
df_NA_nocens <- data.frame(
  time = km_nocens$time,
  H = H_NA_nocens,
  lower = H_NA_nocens - 1.96 * se_H_NA_nocens,
  upper = H_NA_nocens + 1.96 * se_H_NA_nocens
)

# Breslow
km_bres_nocens <- survfit(Surv(time, rep(1,length(time))) ~ 1, type = "kaplan-meier")
df_Bres_nocens <- data.frame(
  time = km_bres_nocens$time,
  H = -log(km_bres_nocens$surv),
  lower = -log(km_bres_nocens$upper),
  upper = -log(km_bres_nocens$lower)
)

# Affichage du graphique
fig <- plot_ly() %>%
  
  # Avec censure
  # Nelson-Aalen
  add_ribbons(data = df_NA_cens, x = ~time, ymin = ~lower, ymax = ~upper,
              fillcolor = "rgba(255,0,0,0.2)", line = list(color = 'transparent'),
              legendgroup = "NA_cens", showlegend = FALSE) %>%
  add_lines(data = df_NA_cens, x = ~time, y = ~H,
            line = list(color = "red", width = 2),
            name = "Nelson-Aalen avec censure",
            legendgroup = "NA_cens", showlegend = TRUE) %>%
  
  # Breslow
  add_ribbons(data = df_Bres_cens, x = ~time, ymin = ~lower, ymax = ~upper,
              fillcolor = "rgba(0,0,255,0.2)", line = list(color = 'transparent'),
              legendgroup = "Bres_cens", showlegend = FALSE) %>%
  add_lines(data = df_Bres_cens, x = ~time, y = ~H,
            line = list(color = "blue", width = 2),
            name = "Breslow avec censure",
            legendgroup = "Bres_cens", showlegend = TRUE) %>%
  
  # Sans censure
  # Nelson-Aalen
  add_ribbons(data = df_NA_nocens, x = ~time, ymin = ~lower, ymax = ~upper,
              fillcolor = "rgba(255,165,0,0.2)", line = list(color = 'transparent'),
              legendgroup = "NA_nocens", showlegend = FALSE) %>%
  add_lines(data = df_NA_nocens, x = ~time, y = ~H,
            line = list(color = "orange", width = 2),
            name = "Nelson-Aalen sans censure",
            legendgroup = "NA_nocens", showlegend = TRUE) %>%
  
  # Breslow
  add_ribbons(data = df_Bres_nocens, x = ~time[-length(df_Bres_nocens$time)],
              ymin = ~lower[-length(df_Bres_nocens$lower)],
              ymax = ~upper[-length(df_Bres_nocens$upper)],
              fillcolor = "rgba(0,128,0,0.2)", line = list(color = 'transparent'),
              legendgroup = "Bres_nocens", showlegend = FALSE) %>%
  add_lines(data = df_Bres_nocens, x = ~time[-length(df_Bres_nocens$time)],
            y = ~H[-length(df_Bres_nocens$H)],
            line = list(color = "green", width = 2),
            name = "Breslow sans censure",
            legendgroup = "Bres_nocens", showlegend = TRUE) %>%
  
  # Layout
  layout(
    title = "Risque cumulatif : avec et sans censure (Nelson-Aalen & Breslow)",
    xaxis = list(title = "Durée du mariage (années)"),
    yaxis = list(title = "Risque cumulatif H(t)"),
    legend = list(orientation = "v", x = 0.1, y = 1)
  )

div(style="text-align:center;", fig)

Les courbes de risque cumulatif sans censure montrent le cumul des divorces en considérant tous les couples comme observés jusqu’à la fin. Pour Nelson-Aalen, le risque cumulatif commence à environ 0,18 après un an, atteint 0,54 après cinq ans et dépasse 5 après quarante ans. Cette courbe reflète le risque brut et tend à surestimer le risque réel, car elle ne tient pas compte des couples censurés. Pour Breslow sans censure, le risque cumulatif est légèrement plus élevé à chaque instant et dépasse 4 après quarante ans, reflétant également une estimation brute du risque de divorce. La courbe ne se termine pas toujours de manière définie à la fin, ce qui est normal puisque le calcul implique le logarithme de la survie et log(0) n’est pas défini.

Les courbes avec censure corrigent le calcul en tenant compte des couples dont on ne connaît pas la fin du mariage. Pour Nelson-Aalen avec censure, le risque cumulatif augmente plus lentement, commençant à environ 0,068 après un an et atteignant environ 2,1 après quarante ans. Les intervalles de confiance s’élargissent progressivement avec le temps car moins de couples restent à risque. Pour Breslow avec censure, le risque cumulatif commence à environ 0,07 après un an et atteint 2,29 après quarante ans.

On voit notamment pour chacune des courbes une augmentation du risque pour la période de 40 ans de durée de marriage. Cela peut s’expliquer de par une plus grande fréquence à 40 ans de durée de marriage parmi les années de 30 ans à 40 ans de durée de marriage selon la distribution de la variable marriage_duration_years.

En résumé, les estimateurs sans censure surestiment le risque cumulatif de divorce, tandis que les estimateurs avec censure donnent une estimation plus réaliste. Nelson-Aalen et Breslow donnent des courbes similaires, Breslow étant légèrement plus lisse et pouvant présenter des valeurs non définies si la survie estimée devient nulle.

3.5.2 - Graphique de la fonction du risque instanné (Nelson-Aalen)

# Données
time <- data.brute$marriage_duration_years
event <- data.brute$divorced

# Hazard instantané Nelson-Aalen avec censure
fit_cens <- survfit(Surv(time, event) ~ 1)
h_cens <- fit_cens$n.event / fit_cens$n.risk
se_h_cens <- sqrt(fit_cens$n.event / (fit_cens$n.risk^2))
df_cens <- data.frame(
  time = fit_cens$time,
  hazard = h_cens,
  lower = h_cens - 1.96 * se_h_cens,
  upper = h_cens + 1.96 * se_h_cens
)

# Hazard instantané Nelson-Aalen sans censure
fit_nocens <- survfit(Surv(time, rep(1, length(time))) ~ 1)
h_nocens <- fit_nocens$n.event / fit_nocens$n.risk
se_h_nocens <- sqrt(fit_nocens$n.event / (fit_nocens$n.risk^2))
df_nocens <- data.frame(
  time = fit_nocens$time,
  hazard = h_nocens,
  lower = h_nocens - 1.96 * se_h_nocens,
  upper = h_nocens + 1.96 * se_h_nocens
)

# Plotly avec IC et legendgroup
fig <- plot_ly() %>%
  # Ribbon Sans censure
  add_ribbons(data = df_nocens, x = ~time, ymin = ~lower, ymax = ~upper,
              fillcolor = 'rgba(0,0,255,0.2)', line = list(color = 'transparent'),
              legendgroup = "nocens", showlegend = FALSE) %>%
  # Courbe centrale Sans censure
  add_trace(data = df_nocens, x = ~time, y = ~hazard, type = 'scatter', mode = 'lines',
            line = list(color = 'blue', width = 2), name = "Sans censure",
            legendgroup = "nocens", showlegend = TRUE) %>%
  # Ribbon Avec censure
  add_ribbons(data = df_cens, x = ~time, ymin = ~lower, ymax = ~upper,
              fillcolor = 'rgba(0,128,0,0.2)', line = list(color = 'transparent'),
              legendgroup = "cens", showlegend = FALSE) %>%
  # Courbe centrale Avec censure
  add_trace(data = df_cens, x = ~time, y = ~hazard, type = 'scatter', mode = 'lines',
            line = list(color = 'darkgreen', width = 2), name = "Avec censure",
            legendgroup = "cens", showlegend = TRUE) %>%
  layout(
    title = "Risque instantané Nelson-Aalen avec IC 95%",
    xaxis = list(title = "Durée du mariage (années)"),
    yaxis = list(title = "Taux de risque h(t)"),
    legend = list(orientation = "v", x = 0.1, y = 1)
  )

div(style = "text-align:center;", fig)

Au début du mariage, le risque instantané de divorce est relativement élevé, avec un pic la première année (0,068 avec censure et 0,178 sans censure). Cela traduit une période initiale plus critique où les couples doivent s’adapter à la vie à deux. Ensuite, le risque se stabilise mais on note des petites augmentations autour de 5 à 7 ans, correspondant possiblement aux premières tensions liées à la vie commune et à l’arrivée éventuelle des enfants.

Une autre période où le risque augmente légèrement se situe entre 25 et 32 ans de mariage, avec des valeurs autour de 0,057–0,058 avec censure et 0,133–0,156 sans censure, reflétant des divorces plus tardifs dus à potentiellement à l’accumulation de tensions sur le long terme.

Enfin, vers la fin de l’observation, la fonction de risque montre des valeurs très élevées, notamment à 40 ans (0,455 avec censure et 1 sans censure), mais il s’agit d’un artefact dû au faible nombre de couples restant à risque et aux calculs de la fonction, et non d’un risque réellement plus élevé dans la population. On observe donc des périodes à risque plus marqué au début, quelques fluctuations intermédiaires et un pic final artificiel.

4 Loi paramétrique : Test de conformité

On cherche à présent à trouver une densité standard permettant d’approcher au mieux la fonction de survie. Pour chacune des distributions standards, on affichera quatre graphiques permettant un diagnostic visuel.

4.1 Loi Gamma

Dans cette section, nous étudions la pertinence de la loi Gamma pour modéliser la durée du mariage (marriage_duration_years), considérée comme une variable de survie positive représentant le temps jusqu’au divorce.

Une variable aléatoire \(T\) suit une loi Gamma de paramètres \(\alpha > 0\) (forme) et \(\theta > 0\) (échelle), notée \(T \sim \Gamma(\alpha, \theta)\), si sa densité de probabilité est donnée par :

\[ f(t) = \frac{1}{\Gamma(\alpha)\theta^\alpha} t^{\alpha-1} e^{-t/\theta}, \quad t > 0 \]

où \(\Gamma(\alpha)\) désigne la fonction Gamma.

Les principaux moments de la loi Gamma sont :

• Espérance : \[ \mathbb{E}[T] = \alpha \theta \]
• Variance : \[ \text{Var}(T) = \alpha \theta^2 \]

---

4.1.1 Interprétation de la loi Gamma en analyse de survie

La loi Gamma est particulièrement adaptée à la modélisation de temps strictement positifs, comme la durée d’un mariage.
Son paramètre de forme \(\alpha\) permet de représenter différentes dynamiques temporelles du risque :

• \(\alpha < 1\) : risque élevé au début puis décroissant (sélection des individus les plus fragiles)
• \(\alpha = 1\) : cas particulier correspondant à une loi exponentielle (risque constant)
• \(\alpha > 1\) : risque croissant dans le temps (accumulation de facteurs de risque)

Contrairement à la loi de Weibull, la loi Gamma ne possède pas une expression simple de la fonction de risque. Elle est donc principalement utilisée pour ajuster correctement la distribution des durées, plutôt que pour interpréter directement l’évolution instantanée du risque. Son intérêt réside dans sa flexibilité et sa capacité à représenter des distributions asymétriques observées empiriquement.

---

4.1.2 Ajustement de la loi Gamma sur la durée du mariage

Dans un premier temps, nous ajustons une loi Gamma à la variable marriage_duration_years afin d’évaluer empiriquement l’adéquation de cette loi aux données observées. Cette étape est purement exploratoire et ne tient pas encore compte de la censure.

# Extraction de la variable de survie
survival_time <- na.omit(data.brute$marriage_duration_years)

# Ajustement de la loi Gamma
fit_gamma <- fitdist(survival_time, "gamma")
summary(fit_gamma)

## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters : 
##       estimate Std. Error
## shape   1.1217    0.01995
## rate    0.1223    0.00272
## Loglikelihood:  -16060   AIC:  32124   BIC:  32137 
## Correlation matrix:
##        shape   rate
## shape 1.0000 0.7997
## rate  0.7997 1.0000

# Résumé des paramètres estimés
plot(fit_gamma,demp=TRUE,histo = TRUE)

4.2 Loi Weibull

Une variable aléatoire \(T\) suit une loi de Weibull de paramètres \(\kappa > 0\) (forme) et \(\lambda > 0\) (échelle), notée \(T \sim \text{Weibull}(\kappa, \lambda)\), si sa densité de probabilité est donnée par :

\[ f(t) = \frac{\kappa}{\lambda} \left(\frac{t}{\lambda}\right)^{\kappa-1} e^{-(t/\lambda)^\kappa}, \quad t > 0 \]

Les principaux moments de la loi de Weibull sont :

• Espérance : \[ \mathbb{E}[T] = \lambda \, \Gamma\left(1 + \frac{1}{\kappa}\right) \]

• Variance : \[ \text{Var}(T) = \lambda^2 \left[ \Gamma\left(1 + \frac{2}{\kappa}\right) - \left(\Gamma\left(1 + \frac{1}{\kappa}\right)\right)^2 \right] \]

où \(\Gamma(\cdot)\) est la fonction Gamma.

---

4.2.1 Interprétation de la loi de Weibull en analyse de survie

La loi de Weibull est très utilisée en analyse de survie car elle permet de représenter différentes dynamiques du risque à travers le paramètre de forme \(\kappa\) :

• \(\kappa < 1\) : risque décroissant dans le temps (événements précoces plus probables)
  
• \(\kappa = 1\) : cas particulier correspondant à une loi exponentielle (risque constant)
  
• \(\kappa > 1\) : risque croissant dans le temps (accumulation des facteurs de risque)

Le paramètre d’échelle \(\lambda\) ajuste la dispersion et le niveau général des durées.
Contrairement à la loi Gamma, la loi de Weibull possède une fonction de risque monotone et est donc particulièrement intéressante pour analyser la dynamique du risque de divorce dans le temps.

---

4.2.2 Ajustement de la loi de Weibull sur la durée du mariage

# Extraction de la variable de survie
survival_time <- na.omit(data.brute$marriage_duration_years)

# Ajustement de la loi de Weibull
fit_weibull <- fitdist(survival_time, "weibull")
summary(fit_weibull)

## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters : 
##       estimate Std. Error
## shape    1.051    0.01151
## scale    9.363    0.13318
## Loglikelihood:  -16070   AIC:  32144   BIC:  32157 
## Correlation matrix:
##        shape  scale
## shape 1.0000 0.3252
## scale 0.3252 1.0000

# Affichage du résumé et des graphiques de diagnostic
plot(fit_weibull, demp=TRUE, histo=TRUE)

4.3 Loi Bêta

Une variable aléatoire \(U\) suit une loi Bêta de paramètres \(\alpha > 0\) et \(\beta > 0\), notée
\(U \sim \text{Beta}(\alpha, \beta)\), si sa densité de probabilité est donnée par :

\[ f(u) = \frac{1}{B(\alpha, \beta)} u^{\alpha-1} (1-u)^{\beta-1}, \quad 0 < u < 1 \]

où \(B(\alpha, \beta)\) est la fonction bêta, définie par :

\[ B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \]

Les principaux moments de la loi Bêta sont :

• Espérance : \[ \mathbb{E}[U] = \frac{\alpha}{\alpha + \beta} \]

• Variance : \[ \text{Var}(U) = \frac{\alpha \beta}{(\alpha + \beta)^2(\alpha + \beta + 1)} \]

---

4.3.1 Interprétation de la loi Bêta en analyse de survie

La loi Bêta n’est pas une loi de survie classique, car elle est définie sur un intervalle borné \([0,1]\).
Elle peut néanmoins être utilisée après normalisation du temps, pour modéliser la position relative d’un événement dans une durée maximale observée.

Les paramètres ont une interprétation intuitive :

• \(\alpha > \beta\) : concentration des événements vers la fin de la durée
• \(\alpha < \beta\) : événements précoces plus fréquents
• \(\alpha = \beta\) : distribution symétrique

La loi Bêta est donc utile pour analyser la structure relative des durées, mais elle ne permet pas d’interpréter directement le risque instantané, contrairement aux lois de Weibull ou exponentielle.

---

4.3.2 Ajustement de la loi Bêta sur la durée du mariage normalisée

Avant l’ajustement, la variable marriage_duration_years est normalisée sur l’intervalle \([0,1]\) en la divisant par la durée maximale observée.

# Extraction de la variable de survie
survival_time <- na.omit(data.brute$marriage_duration_years)

# Normalisation sur [0,1]
survival_time_scaled <- survival_time / max(survival_time)

# Ajustement de la loi Bêta
fit_beta <- fitdist(survival_time_scaled, "beta")

# Affichage des graphiques de diagnostic
summary(fit_beta)

## Fitting of the distribution ' beta ' by maximum likelihood 
## Parameters : 
##        estimate
## shape1   0.5989
## shape2   2.0131
## Loglikelihood:  -1.074e+13   AIC:  2.147e+13   BIC:  2.147e+13

plot(fit_beta, demp=TRUE, histo=TRUE)

4.4 Loi exponentielle

Une variable aléatoire \(T\) suit une loi exponentielle de paramètre \(\lambda > 0\), notée
\(T \sim \text{Exp}(\lambda)\), si sa densité de probabilité est donnée par :

\[ f(t) = \lambda e^{-\lambda t}, \quad t > 0 \]

Les principaux moments de la loi exponentielle sont :

• Espérance : \[ \mathbb{E}[T] = \frac{1}{\lambda} \]

• Variance : \[ \text{Var}(T) = \frac{1}{\lambda^2} \]

---

4.4.1 Interprétation de la loi exponentielle en analyse de survie

La loi exponentielle est le modèle de survie le plus simple.
Elle repose sur l’hypothèse clé d’un risque constant dans le temps :

\[ h(t) = \lambda \]

Cela signifie que la probabilité instantanée de divorce est indépendante de la durée déjà écoulée du mariage.

Cette propriété implique également la propriété de perte de mémoire :

\[ P(T > t+s \mid T > t) = P(T > s) \]

En pratique, cette hypothèse est souvent forte dans le contexte du divorce, car elle suppose que : - un mariage récent et un mariage ancien ont le même risque instantané de rupture.

La loi exponentielle est donc principalement utilisée comme : - modèle de référence (baseline), - ou point de comparaison pour des modèles plus flexibles (Weibull, Cox).

---

# Extraction de la variable de survie
survival_time <- na.omit(data.brute$marriage_duration_years)

# Ajustement de la loi exponentielle
fit_exp <- fitdist(survival_time, "exp")
summary(fit_exp)

## Fitting of the distribution ' exp ' by maximum likelihood 
## Parameters : 
##      estimate Std. Error
## rate    0.109   0.001542
## Loglikelihood:  -16080   AIC:  32162   BIC:  32169

# Affichage des graphiques de diagnostic
plot(fit_exp, demp=TRUE, histo=TRUE)

4.5 Loi normale

Une variable aléatoire \(T\) suit une loi normale de paramètres \(\mu \in \mathbb{R}\) (moyenne) et \(\sigma > 0\) (écart-type), notée
\(T \sim \mathcal{N}(\mu, \sigma^2)\), si sa densité de probabilité est donnée par :

\[ f(t) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right), \quad t \in \mathbb{R} \]

Les principaux moments de la loi normale sont :

• Espérance : \[ \mathbb{E}[T] = \mu \]

• Variance : \[ \text{Var}(T) = \sigma^2 \]

---

4.5.1 Interprétation de la loi normale en analyse de survie

La loi normale n’est pas adaptée en théorie à l’analyse de survie, car :

• elle autorise des durées négatives,
• elle ne respecte pas la contrainte de positivité des temps,
• elle ne possède pas de fonction de survie ou de risque interprétable dans ce contexte.

Cependant, elle peut être utilisée à titre exploratoire pour :

• vérifier si la distribution empirique des durées est approximativement symétrique,
• comparer visuellement la forme des données à une distribution classique,
• motiver l’utilisation de transformations (logarithme, par exemple).

Dans le contexte du divorce, une bonne adéquation à une loi normale serait généralement peu plausible, car les durées de mariage sont souvent asymétriques et bornées inférieurement par zéro.

---

Nous ajustons ici une loi normale à la variable marriage_duration_years à des fins purement exploratoires, afin de comparer sa forme à celle des lois de survie paramétriques usuelles.

# Extraction de la variable de survie
survival_time <- na.omit(data.brute$marriage_duration_years)

# Ajustement de la loi normale
fit_norm <- fitdist(survival_time, "norm")
summary(fit_norm)

## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##      estimate Std. Error
## mean    9.171    0.12512
## sd      8.847    0.08847
## Loglikelihood:  -17995   AIC:  35995   BIC:  36008 
## Correlation matrix:
##      mean sd
## mean    1  0
## sd      0  1

# Affichage des graphiques de diagnostic
plot(fit_norm, demp=TRUE, histo=TRUE)

Avec les résulats obtenu précédement nous n’exploreront pas la loi lognormale qui ne semble pas pertinente ici.

4.6 Comparaison des distributions paramétriques par AIC, BIC et test du chi-carré

4.6.1 Rappel sur les critères AIC et BIC

Les critères AIC (Akaike Information Criterion) et BIC (Bayesian Information Criterion) permettent de comparer plusieurs distributions paramétriques ajustées aux données, en prenant en compte :

• la qualité de l’ajustement (vraisemblance),
• la complexité du modèle (nombre de paramètres).

Une valeur plus faible indique un meilleur compromis entre ajustement et parcimonie.

---

4.6.2 Test du chi-carré d’adéquation

Le test du chi-carré d’adéquation permet d’évaluer si une distribution théorique ajustée est compatible avec les données observées.

• Hypothèse nulle \(H_0\) :
  > Les données suivent la distribution théorique considérée.

• Hypothèse alternative \(H_1\) :
  > Les données ne suivent pas cette distribution.

Une p-value faible conduit au rejet de \(H_0\), indiquant une mauvaise adéquation du modèle.
Une p-value élevée suggère que l’on ne peut pas rejeter l’hypothèse d’adéquation.

---

4.6.3 Tableau comparatif des lois testées

Comparaison des lois paramétriques par AIC, BIC et test du chi-carré

Distribution	AIC	BIC	Chi-deux	p-value
Gamma	3.212e+04	3.214e+04	621.6	< 2.2e-16
Weibull	3.214e+04	3.216e+04	500.3	< 2.2e-16
Exponentielle	3.216e+04	3.217e+04	332.1	< 2.2e-16
Normale	3.599e+04	3.601e+04	3136.8	< 2.2e-16
Bêta (normalisée)	2.147e+13	2.147e+13	919.8	< 2.2e-16

5 Modèle Cox

On se place dans le cadre du modèle de Cox. La fonction de risque conditionnelle s’écrit

\[ h(t \mid X) = h_0(t)\exp(\beta^T X), \]

où :

• \(X = (X_1, X_2, \dots, X_p)^T\) est le vecteur de covariables de dimension \(p \times 1\),
• \(\beta = (\beta_1, \beta_2, \dots, \beta_p)^T\) est le vecteur de coefficients de régression.
• \(h_0(t)\) est le risque instantané d’évènement quand toutes les covariables valent zéros.

Considérons :

• \(D\) le nombre d’événements divorce observés parmi les \(n\) sujets de l’étude,
• \(T_1 < T_2 < \dots < T_D\) les temps distincts de divorce,
• \(R(T_i)\) l’ensemble des individus encore à risque de divorce juste avant le temps \(T_i\) (toujours mariés à \(T_i^-\)).

Si l’on compare le risque instantané d’un individu ayant \(X_j = 1\) à celui d’un individu de référence dont toutes les covariables valent zéro (\(X = 0\)), on a :

\[ \frac{h\bigl(t \,\big|\, X_j = 1, \, X_i = 0 \ \forall i \neq j\bigr)}{h_0(t \mid X = 0)} = \exp(\beta_j) = \text{HR}_j \]

En pratique, le HR indique combien de fois le risque instantané est multiplié pour une augmentation d’une unité de la covariable \(X_j\), comparé à un individu de référence avec toutes les autres covariables à zéro.

5.1 Tests d’hypothèses globaux

À partir de la log-vraisemblance partielle, on peut estimer le vecteur de paramètres β de dimension \(p \times 1\) :

\[ L(\beta) = \log L_{\text{Cox}}(\beta) = \sum_{i=1}^{D} \left[ \beta^T X_{(i)} - \log \left( \sum_{j \in R(T_i)} \exp(\beta^T X_j) \right) \right], \]

La fonction score \(U(\beta)\) est le vecteur des dérivées premières de la log-vraisemblance :

\[ U(\beta) = \frac{\partial L(\beta)}{\partial \beta} = \sum_{i=1}^{D} \left[ X_{(i)} - \frac{\sum_{j \in R(T_i)} X_j \exp(\beta^T X_j)}{\sum_{j \in R(T_i)} \exp(\beta^T X_j)} \right] = \begin{pmatrix} \frac{\partial L(\beta)}{\partial \beta_1} \\[1mm] \vdots \\[1mm] \frac{\partial L(\beta)}{\partial \beta_p} \end{pmatrix}. \]

Un estimateur consistant de la matrice de variance-covariance de \(\hat{\beta}\) est donné par l’inverse de la matrice d’information de Fisher :

\[ \mathrm{Var}(\hat{\beta}) = \hat{I}(\hat{\beta})^{-1}, \quad \text{avec } [I(\beta)]_{i,j} = - \frac{\partial^2 L(\beta)}{\partial \beta_i \partial \beta_j}. \]

On souhaite tester :

\[ \begin{cases} H_0: \beta_1 = \dots = \beta_p = 0 \\[1mm] H_1: \exists j \in \{1,\dots,p\} \text{ tel que } \beta_j \neq 0 \end{cases} \]

Trois statistiques sont classiquement utilisées :

1. Test du rapport de vraisemblance (Likelihood Ratio Test) :

\[ \chi^2_{\text{LRT}} = 2 \left[ \log L_{\text{Cox}}(\hat{\beta}) - \log L_{\text{Cox}}(\beta_0) \right] \sim \chi^2_p \text{ sous } H_0 \]

2. Test de Wald :

\[ \chi^2_W = (\hat{\beta} - \beta_0)^T I(\hat{\beta}) (\hat{\beta} - \beta_0) \sim \chi^2_p \text{ sous } H_0 \]

3. Test de Score (ou logrank) :

\[ \chi^2_S = U(\beta_0)^T I(\beta_0)^{-1} U(\beta_0) \sim \chi^2_p \text{ sous } H_0 \]

Toutes ces statistiques suivent une loi du chi-deux à \(p\) degrés de liberté sous l’hypothèse nulle. Elles permettent de déterminer si au moins une covariable a un effet significatif sur le risque de l’événement. De plus \(\beta_0\) = 0 dans cette section.

# Données
time <- data.brute$marriage_duration_years
event <- data.brute$divorced

# Modèle de Cox
cox_all <- coxph(
  Surv(time, event) ~ . - divorced - marriage_duration_years,
  data = data.brute
)

model_cox <- summary(cox_all)

# Résultats des tests globaux
cat(
  "Likelihood ratio test =", model_cox$logtest[1], "on", model_cox$logtest[2], "df, p=", signif(model_cox$logtest[3], 3), "\n",
  "Wald test            =", model_cox$waldtest[1], "on", model_cox$waldtest[2], "df, p=", signif(model_cox$waldtest[3], 3), "\n",
  "Score (logrank) test =", model_cox$sctest[1], "on", model_cox$sctest[2], "df, p=", signif(model_cox$sctest[3], 3), "\n"
)

## Likelihood ratio test = 117.1 on 29 df, p= 1.54e-12 
##  Wald test            = 120.4 on 29 df, p= 4.21e-13 
##  Score (logrank) test = 120.7 on 29 df, p= 3.71e-13

Le modèle de Cox ajusté pour les covariables étudiées est globalement significatif. Les trois tests globaux Likelihood ratio, Wald et Score (logrank) donnent des statistiques élevées avec des p-values très faibles (respectivement 1.54e-12, 4.21e-13 et 3.71e-13). Sous l’hypothèse nulle, ces tests évaluent si tous les coefficients des covariables sont nuls, c’est-à-dire si aucune covariable n’a d’effet sur le risque de divorce. Le rejet de cette hypothèse nulle indique qu’au moins une covariable a un effet significatif sur le risque de divorce, justifiant l’examen des coefficients individuels pour identifier quelles covariables contribuent de manière significative au modèle.

5.2 Tests par covariable

La colonne Pr(>|z|) correspond au test de Wald pour chaque coefficient individuel.

Chaque covariable a sa propre hypothèse nulle :

\[ H_0 : \beta_j = 0 \quad \text{vs} \quad H_1 : \beta_j \neq 0 \]

La statistique utilisée est la statistique de Wald par covariable :

\[ z_j = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)}, \quad \text{p-value} = 2 \cdot \Pr(|Z| > |z_j|), \]

où \(\hat{\beta}_j\) est l’estimation du coefficient et \(\text{SE}(\hat{\beta}_j)\) son erreur-type tel que \(\text{SE}(\hat{\beta}_j) = \sqrt{\text{Var}(\hat{\beta}_j)} = \sqrt{\left[ \hat{I}(\hat{\beta})^{-1} \right]_{jj}},\) .

# Données
time <- data.brute$marriage_duration_years
event <- data.brute$divorced

cox_all <- coxph(
  Surv(time, event) ~ . - divorced - marriage_duration_years,
  data = data.brute
)

model_cox <- summary(cox_all)

model_cox$coefficients

##                                              coef exp(coef)  se(coef)        z
## age_at_marriage                         5.486e-03    1.0055 4.605e-03  1.19132
## num_children                            4.575e-03    1.0046 1.829e-02  0.25018
## education_levelHigh School              5.819e-02    1.0599 5.409e-02  1.07579
## education_levelMaster                   7.091e-02    1.0735 6.151e-02  1.15270
## education_levelNo Formal Education     -1.521e-01    0.8589 1.178e-01 -1.29154
## education_levelPhD                      1.685e-01    1.1835 1.058e-01  1.59210
## employment_statusHomemaker              5.806e-02    1.0598 6.586e-02  0.88164
## employment_statusPart-time              1.587e-02    1.0160 6.021e-02  0.26349
## employment_statusUnemployed             4.198e-02    1.0429 6.619e-02  0.63423
## combined_income                         1.160e-06    1.0000 1.137e-06  1.01971
## religious_compatibilityNot Religious   -1.642e-01    0.8486 7.232e-02 -2.27006
## religious_compatibilitySame Religion    4.055e-03    1.0041 5.872e-02  0.06905
## cultural_background_match               1.955e-02    1.0197 5.217e-02  0.37476
## communication_score                    -4.817e-02    0.9530 1.130e-02 -4.26250
## conflict_frequency                      6.002e-03    1.0060 1.576e-02  0.38073
## conflict_resolution_styleAvoidant      -9.865e-02    0.9061 6.823e-02 -1.44582
## conflict_resolution_styleCollaborative -5.388e-02    0.9475 6.157e-02 -0.87503
## conflict_resolution_stylePassive       -2.855e-03    0.9971 7.616e-02 -0.03749
## financial_stress_level                  2.739e-02    1.0278 9.770e-03  2.80366
## mental_health_issues                    1.602e-01    1.1738 5.424e-02  2.95406
## infidelity_occurred                     2.389e-01    1.2698 5.863e-02  4.07435
## counseling_attended                     5.089e-02    1.0522 5.247e-02  0.96983
## social_support                         -1.708e-02    0.9831 1.131e-02 -1.51070
## shared_hobbies_count                   -2.011e-02    0.9801 1.342e-02 -1.49929
## marriage_typeLove                       2.732e-02    1.0277 5.304e-02  0.51511
## marriage_typeOther                     -1.484e-01    0.8621 1.134e-01 -1.30920
## pre_marital_cohabitation               -5.820e-02    0.9435 4.583e-02 -1.26972
## domestic_violence_history               3.665e-01    1.4426 8.761e-02  4.18290
## trust_score                            -5.273e-02    0.9486 1.145e-02 -4.60496
##                                         Pr(>|z|)
## age_at_marriage                        2.335e-01
## num_children                           8.024e-01
## education_levelHigh School             2.820e-01
## education_levelMaster                  2.490e-01
## education_levelNo Formal Education     1.965e-01
## education_levelPhD                     1.114e-01
## employment_statusHomemaker             3.780e-01
## employment_statusPart-time             7.922e-01
## employment_statusUnemployed            5.259e-01
## combined_income                        3.079e-01
## religious_compatibilityNot Religious   2.320e-02
## religious_compatibilitySame Religion   9.450e-01
## cultural_background_match              7.078e-01
## communication_score                    2.022e-05
## conflict_frequency                     7.034e-01
## conflict_resolution_styleAvoidant      1.482e-01
## conflict_resolution_styleCollaborative 3.816e-01
## conflict_resolution_stylePassive       9.701e-01
## financial_stress_level                 5.053e-03
## mental_health_issues                   3.136e-03
## infidelity_occurred                    4.614e-05
## counseling_attended                    3.321e-01
## social_support                         1.309e-01
## shared_hobbies_count                   1.338e-01
## marriage_typeLove                      6.065e-01
## marriage_typeOther                     1.905e-01
## pre_marital_cohabitation               2.042e-01
## domestic_violence_history              2.878e-05
## trust_score                            4.126e-06

En regardant donc les coefficients, on observe de nombreuses variables ayant un effet significatif sur le modèle.

La compatibilité religieuse apparaît comme un facteur protecteur : les couples où la personne n’est pas religieuse ont un hazard ratio de 0.849, ce qui correspond à une réduction du risque de divorce d’environ 15 % par rapport à la catégorie de référence. Cet effet est statistiquement significatif (z = -2.27, p = 0.0232), indiquant que l’absence de religion réduit légèrement mais de manière significative le risque de rupture.

La qualité de la communication entre partenaires a un effet protecteur important. Chaque augmentation d’une unité du score de communication est associée à un hazard ratio de 0.953, soit une diminution du risque de divorce d’environ 4,7 %, avec une très forte significativité statistique (z = -4.26, p = 2.0e-05). Cela indique que de meilleures compétences en communication réduisent le risque de divorce.

Le stress financier augmente le risque de divorce : chaque unité supplémentaire dans le niveau de stress financier est associée à un hazard ratio de 1.028, soit une augmentation du risque de divorce d’environ 2,8 %, avec z = 2.80 et p = 0.0051.

Les problèmes de santé mentale accroissent également le risque de divorce. Les couples où un partenaire présente des problèmes de santé mentale ont un hazard ratio de 1.174, ce qui correspond à une augmentation du risque de divorce d’environ 17,4 %, z = 2.95, p = 0.0031.

L’infidélité est un facteur de risque majeur : les couples où l’infidélité est survenue ont un hazard ratio de 1.270, soit une augmentation du risque de divorce d’environ 27 %, z = 4.07, p = 4.6e-05.

L’historique de violence domestique est le facteur le plus impactant observé : les couples ayant un tel historique présentent un hazard ratio de 1.443, correspondant à une augmentation du risque de divorce d’environ 44,3 %, avec z = 4.18 et p = 2.9e-05.

Enfin, le niveau de confiance est fortement protecteur : chaque unité supplémentaire dans le score de confiance est associée à un hazard ratio de 0.949, soit une réduction du risque de divorce d’environ 5,1 %, avec une significativité très élevée (z = -4.60, p = 4.1e-06).

Ces résultats montrent que les facteurs relationnels (communication, confiance, infidélité, violence domestique), psychologiques (problèmes de santé mentale) et financiers (stress financier) ont un impact significatif sur la stabilité du mariage. Les variables avec HR < 1 réduisent le risque de divorce, tandis que celles avec HR > 1 l’augmentent.

5.3 Synthèse statistique du modèle Cox

# Création du dataframe synthétique
df_significatives <- data.frame(
  Variable = c(
    "Religious Compatibility (Not Religious)",
    "Communication Score",
    "Financial Stress Level",
    "Mental Health Issues",
    "Infidelity Occurred",
    "Domestic Violence History",
    "Trust Score"
  ),
  HR = c(0.849, 0.953, 1.028, 1.174, 1.270, 1.443, 0.949),
  Pourcentage_Risque = c(
    "-15%", "-4.7%", "+2.8%", "+17.4%", "+27%", "+44.3%", "-5.1%"
  ),
  Interpretation = c(
    "Non religieuse → réduction du risque de divorce",
    "Meilleure communication → réduction du risque",
    "Stress financier ↑ → risque de divorce ↑",
    "Problèmes de santé mentale ↑ → risque ↑",
    "Infidélité → risque de divorce ↑",
    "Historique de violence → risque ↑",
    "Confiance ↑ → réduction du risque de divorce"
  )
)

# Coloration dynamique de la colonne Pourcentage_Risque
df_significatives$Pourcentage_Risque <- sapply(df_significatives$Pourcentage_Risque, function(x) {
  value <- as.numeric(gsub("[+%]", "", x))  # extraire valeur numérique
  color <- ifelse(value > 0, "red", "green")
  cell_spec(x, color = color, bold = TRUE)
})

# Création du tableau
kable(df_significatives, "html", escape = FALSE,
      caption = "📊 Covariables significatives du modèle de Cox") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  column_spec(1, bold = TRUE, color = "black", background = c("#D6EAF8", "#F9E79F", "#FDEDEC", "#D6EAF8", "#F9E79F", "#FDEDEC", "#D6EAF8")) %>%
  column_spec(2, width = "3cm") %>%
  column_spec(3, width = "3cm") %>%
  column_spec(4, width = "8cm") %>%
  row_spec(0, bold = TRUE, color = "white", background = "#2C3E50")

📊 Covariables significatives du modèle de Cox

Variable	HR	Pourcentage_Risque	Interpretation
Religious Compatibility (Not Religious)	0.849	-15%	Non religieuse → réduction du risque de divorce
Communication Score	0.953	-4.7%	Meilleure communication → réduction du risque
Financial Stress Level	1.028	+2.8%	Stress financier ↑ → risque de divorce ↑
Mental Health Issues	1.174	+17.4%	Problèmes de santé mentale ↑ → risque ↑
Infidelity Occurred	1.270	+27%	Infidélité → risque de divorce ↑
Domestic Violence History	1.443	+44.3%	Historique de violence → risque ↑
Trust Score	0.949	-5.1%	Confiance ↑ → réduction du risque de divorce

6 Conclusion

À travers notre analyse de survie, nous avons pu observer diverses applications des fonctions de survie dans l’étude du divorce. La comparaison entre les données censurées et non censurées met en évidence des différences notables dans l’estimation de la survie. Certaines covariables semblent particulièrement influentes, comme les problèmes de santé mentale, qui montrent un impact significatif sur la durée de vie d’un mariage. En revanche, le niveau d’éducation ne semble pas produire de différence significative.

L’utilisation des données non censurées pose des limites pour les estimateurs de Breslow de la fonction de risque, ce qui suggère de privilégier l’estimateur de Nelson-Aalen pour ces données. Néanmoins, les deux méthodes donnent des résultats globalement similaires, à une asymptote près.

Notre analyse révèle que certains facteurs influencent fortement le risque de divorce, notamment durant la première année de mariage, avec une augmentation notable du risque entre 25 et 30 ans de durée de mariage. Selon les critères AIC et BIC, la loi gamma représente le meilleur ajustement, avec des paramètres estimés à \(\alpha\) = 1.1217 et \(\theta\) = 0.1223.

Enfin, le modèle de Cox met en évidence l’influence de variables clés sur la survie d’un mariage : une bonne communication réduit le risque de divorce, tandis que la présence de violence domestique l’augmente considérablement.

7 Référence

• Saint Pierre, P. (2021). Introduction à l’analyse des durées de survie. Institut de Mathématiques de Toulouse, Université Paul Sabatier.